# Loading required packages
library(glmmTMB); library(lme4); library(parameters); library(Matrix); library(methods); library(mgcv); library(nlme); library(numDeriv); library(stargazer); library(tidyverse); library(TMB)
# Loading the data
load("county_level_request_data.RData")

#############
# TABLE A.3 #
#############
# Making a table of means and standard deviations for county-level request (non-standardized) variables
# Variables to use
cols <- c('log_sum', "dem", "on_appropriations", "meddist", "log_countypop", "other_sen_requested", "other_sen_sameparty", "female", 
          "core_county","swing_county","seniority","party_leader", "freshman", "pct_urban", "median_household_income", "capital")

stargazer(as.data.frame(senators_appropriations[,cols]), covariate.labels=c("Log(County Appropriations Requests + 1)", "Democrat (Majority Party Member)", 
                                                                            "Member of Appropriations Committee","Distance from Floor Median",
                                                                            "Log County Population","Other Senator Requested Funds to County","Other Senator is Same Party",
                                                                            "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                                                                            "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                                                                            "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          title="Summary Statistics of Dataset Analyzing County-Level Outcomes", digits=2, summary.stat= c("mean", "sd", "min","max"),
          notes='\\parbox[t]{1\\textwidth}{\\footnotesize \\textit{Note}: Table presents summary statistics for the (non-standardized) variables used in our county-level analysis.}')

#############
# TABLE C.2 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# FY 2022
# Running the zero inflated model
zim_2_fy2022 <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital,
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[senators_appropriations$year==2021,], family = gaussian)
# Model summary and standard errors
summary_zim_fy2022 <- summary(zim_2_fy2022)
ses_fy2022 <-standard_error(zim_2_fy2022)

# Save coefficients from the model
cond_fy2022 <- summary_zim_fy2022$coefficients$cond
zi_fy2022 <- summary_zim_fy2022$coefficients$zi
coefs_fy2022 <- c(cond_fy2022[,1])
coefs2_fy2022 <- c(zi_fy2022[,1])

# Save standard errors from the model
ses_cond_fy2022 <- c(ses_fy2022$SE[ses_fy2022$Component=="conditional"])
names(ses_cond_fy2022) <- ses_fy2022$Parameter[ses_fy2022$Component=="conditional"]
ses_zi_fy2022 <- c(ses_fy2022$SE[ses_fy2022$Component=="zero_inflated"])
names(ses_zi_fy2022) <- ses_fy2022$Parameter[ses_fy2022$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_fy2022 <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital,
                data = senators_appropriations[senators_appropriations$year==2021,])

basic_reg2_fy2022 <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[senators_appropriations$year==2021 & complete.cases(senators_appropriations),])

# FY 2023
# Running the zero inflated model
zim_fy_2023 <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                      core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                      freshman + pct_urban_scaled + median_household_income_scaled + capital,
                    zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                    data = senators_appropriations[senators_appropriations$year==2022,], family = gaussian)
# Model summary and standard errors
summary_zim_fy_2023 <- summary(zim_fy_2023)
ses_fy_2023 <-standard_error(zim_fy_2023)

# Save coefficients from the model
cond_fy_2023 <- summary_zim_fy_2023$coefficients$cond
zi_fy_2023 <- summary_zim_fy_2023$coefficients$zi
coefs_fy_2023 <- c(cond_fy_2023[,1])
coefs2_fy_2023 <- c(zi_fy_2023[,1])

# Save standard errors from the model
ses_cond_fy_2023 <- c(ses_fy_2023$SE[ses_fy_2023$Component=="conditional"])
names(ses_cond_fy_2023) <- ses_fy_2023$Parameter[ses_fy_2023$Component=="conditional"]
ses_zi_fy_2023 <- c(ses_fy_2023$SE[ses_fy_2023$Component=="zero_inflated"])
names(ses_zi_fy_2023) <- ses_fy_2023$Parameter[ses_fy_2023$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_fy_2023 <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                       other_sen_requested*other_sen_sameparty +
                       female + core_county + swing_county + seniority_scaled + party_leader + 
                       freshman + pct_urban_scaled + median_household_income_scaled + capital,
                     data = senators_appropriations[senators_appropriations$year==2022,])

basic_reg2_fy_2023 <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                      data = senators_appropriations[senators_appropriations$year==2022 & complete.cases(senators_appropriations),])


# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
pooled.model.list <- list(basic_reg2_fy2022, basic_reg_fy2022, basic_reg2_fy_2023, basic_reg_fy_2023)
stargazer(pooled.model.list,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage 2022","Second Stage 2022","First Stage 2023","Second Stage 2023"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Other Sen Requested * Other Sen Same Party"),
          notes="\\parbox[t]{\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_sepyears",
          digits=3,
          coef=list(coefs2_fy2022, coefs_fy2022, coefs2_fy_2023, coefs_fy_2023), 
          se=list(ses_zi_fy2022, ses_cond_fy2022, ses_zi_fy_2023, ses_cond_fy_2023),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Subset by Fiscal Year",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

#############
# TABLE C.5 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Pooling both years together
senators_appropriations_oneyear <- senators_appropriations %>%
  group_by(senator,state,GEOID, log_countypop, on_appropriations,female, dem,
           freshman, core_county, swing_county, seniority, party_leader, meddist,
           capital, pct_urban, median_household_income)%>%
  dplyr::summarize(log_sum = log(sum(appropriation_sum)+1),
                   other_sen_requested = max(other_sen_requested),
                   other_sen_sameparty = max(other_sen_sameparty))

# Re-standardize the non-binary variables (now 6284 observations instead of 12568)
# Variables to standardize
senators_appropriations_oneyear_scaled <- senators_appropriations_oneyear[,c("log_countypop", "seniority", "meddist", "pct_urban", "median_household_income")]
# Standardize
senators_appropriations_oneyear_scaled <- as.data.frame(scale(senators_appropriations_oneyear_scaled))
# Rename standardized variables
colnames(senators_appropriations_oneyear_scaled) <- c("log_countypop_scaled", "seniority_scaled", "meddist_scaled", "pct_urban_scaled", "median_household_income_scaled")
# Add the standardized variables to the dataset
senators_appropriations_oneyear <- bind_cols(senators_appropriations_oneyear, senators_appropriations_oneyear_scaled)

# Running the zero inflated model
zim_2_oneyear_county <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                                  core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                                  freshman + pct_urban_scaled + median_household_income_scaled + capital,
                                zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                                data = senators_appropriations_oneyear, family =  gaussian)

# Model summary and standard errors
summary_zim_oneyear_county <- summary(zim_2_oneyear_county)
ses_oneyear_county <- standard_error(zim_2_oneyear_county)

# Save coefficients from the model
cond_oneyear_county <- summary_zim_oneyear_county$coefficients$cond
zi_oneyear_county <- summary_zim_oneyear_county$coefficients$zi
coefs_oneyear_county <- c(cond_oneyear_county[,1])
coefs2_oneyear_county <- c(zi_oneyear_county[,1])

# Save standard errors from the model
ses_cond_oneyear_county <- c(ses_oneyear_county$SE[ses_oneyear_county$Component=="conditional"])
names(ses_cond_oneyear_county) <- ses_oneyear_county$Parameter[ses_oneyear_county$Component=="conditional"]
ses_zi_oneyear_county <- c(ses_oneyear_county$SE[ses_oneyear_county$Component=="zero_inflated"])
names(ses_zi_oneyear_county) <- ses_oneyear_county$Parameter[ses_oneyear_county$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_oneyear_county <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital,
                data = senators_appropriations_oneyear)

basic_reg2_oneyear_county <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations_oneyear[complete.cases(senators_appropriations_oneyear),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_oneyear_county, basic_reg_oneyear_county,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_pooled",
          digits=3,
          coef=list(coefs2_oneyear_county, coefs_oneyear_county), 
          se=list(ses_zi_oneyear_county, ses_cond_oneyear_county),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Pooling Fiscal Years",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

#############
# TABLE C.7 #
#############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_main <- glmmTMB(log_sum ~ log_countypop_scaled  + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                    core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                    freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                  zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                  data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_main <- summary(zim_main)
ses_main <- standard_error(zim_main)

# Save coefficients from the model
cond_main <- summary_zim_main$coefficients$cond
zi_main <- summary_zim_main$coefficients$zi
coefs_main <- c(cond_main[,1])
coefs2_main <- c(zi_main[,1])

# Save standard errors from the model
ses_cond_main <- c(ses_main$SE[ses_main$Component=="conditional"])
names(ses_cond_main) <- ses_main$Parameter[ses_main$Component=="conditional"]
ses_zi_main <- c(ses_main$SE[ses_main$Component=="zero_inflated"])
names(ses_zi_main) <- ses_main$Parameter[ses_main$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_main <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled  +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                data = senators_appropriations)

basic_reg2_main <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_main, basic_reg_main,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county",
          digits=3,
          coef=list(coefs2_main, coefs_main), 
          se=list(ses_zi_main, ses_cond_main),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.10 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Set the seed to run the bootstrap and create empty lists to store results
set.seed(02138)
# 20000 bootstraps
nboot <- 20000
results_cond_coefs  <- list()
results_cond_ses  <- list()
results_zi_coefs  <- list()
results_zi_ses  <- list()

# Creating a unique senator variable
senators_appropriations$unique_senator_name <- paste0(senators_appropriations$senator, senators_appropriations$state)

# Bootstrapping at the senator level
ids <- unique(senators_appropriations$unique_senator_name)

# Running the bootstrap
# Ensure R and package versions match readme
for (k in 1:nboot) {
  # Sample senators to be used in this bootstrap iteration
  senator_sample <- sample(ids,  length(ids),replace=TRUE)
  
  # Take all of the counties that each observation in our bootstrapped sample has
  dataset <- inner_join(tibble(unique_senator_name = senator_sample), senators_appropriations, by = "unique_senator_name")
  
  # Then, run the zero-inflated model with this sample
  zim_2 <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                     core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                     freshman + pct_urban_scaled + median_household_income_scaled + factor(year) + capital,
                   zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                   data = dataset, family =  gaussian)

  summary_zim <- summary(zim_2)
  
  # Storing the results
  results_cond_coefs[[k]]  <- summary_zim$coefficients$cond[,1]
  results_cond_ses[[k]]  <- summary_zim$coefficients$cond[,2]
  results_zi_coefs[[k]]  <- summary_zim$coefficients$zi[,1]
  results_zi_ses[[k]]  <- summary_zim$coefficients$zi[,2]
  
  # Printing progress
  if (k %% 100 == 0) cat(paste(k, "out of", nboot, "samples stored..\n"))
}

# Turn the bootstrap results into our CI estimates
# Omitting models with NA SEs
to.keep <- which(apply(is.na(do.call(rbind.data.frame,results_cond_ses)),1,sum) == 0)

# SD of coefficients to create standard errors for bootstrap
cond_vec_se_alt <- apply(do.call(rbind.data.frame,results_cond_coefs)[to.keep,], 2, sd)
zi_vec_se_alt <- apply(do.call(rbind.data.frame,results_zi_coefs)[to.keep,], 2, sd)

# Variable names for our bootstrapped SEs
names(cond_vec_se_alt) <- c("intercept","log_countypop_scaled", "other_sen_requested", "other_sen_sameparty", "on_appropriations",
                            "female", "dem", "core_county","swing_county","seniority_scaled",
                            "party_leader", "meddist_scaled", "freshman", "pct_urban_scaled", 
                            "median_household_income_scaled", "factor(year)2022","capital","other_sen_requested:other_sen_sameparty")

names(zi_vec_se_alt) <- c("intercept","dem", "on_appropriations", "meddist_scaled", "dem:meddist_scaled")

# Rounding
cond_vec_se_alt <- round(cond_vec_se_alt,3)
zi_vec_se_alt <- round(zi_vec_se_alt,3)

# Running the zero inflated model outside of the loop to get coefficients
zim_2_boot <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary
summary_zim_boot <- summary(zim_2_boot)

# Save coefficients from the model
cond_boot <- summary_zim_boot$coefficients$cond
zi_boot <- summary_zim_boot$coefficients$zi
coefs_boot <- c(cond_boot[,1])
coefs2_boot <- c(zi_boot[,1])

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model and the bootstrap above
basic_reg_boot <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + factor(year) + capital,
                data = senators_appropriations)

basic_reg2_boot <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# se takes the bootstrapped SEs
stargazer(basic_reg2_boot, basic_reg_boot,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from DW-NOMINATE Median",
                              "Distance from DW-NOMINATE Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "Fiscal Year 2023","County Includes Capital City", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. Standard errors come from 19491 bootstraps conducted at the senator level. 
                        We are forced to drop 509 of our initial 20000 bootstrapped iterations due to insufficient coverage for some of our explanatory variables. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_boot",
          digits=3,
          coef=list(coefs2_boot, coefs_boot), 
          se=list(zi_vec_se_alt, cond_vec_se_alt),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level with Bootstrapped Standard Errors",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.11 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Set the seed to run the bootstrap and create empty lists to store results
set.seed(02138)
# 20000 bootstraps
nboot <- 20000
results_cond_coefs  <- list()
results_cond_ses  <- list()
results_zi_coefs  <- list()
results_zi_ses  <- list()

# Creating a unique senator variable
senators_appropriations$unique_senator_name <- paste0(senators_appropriations$senator, senators_appropriations$state)

# Bootstrapping at the county level
ids <- unique(senators_appropriations$GEOID)

# Running the bootstrap
# Ensure R and package versions match readme
for (k in 1:nboot) {
  # Sample counties to be used in this bootstrap iteration
  county_sample <- sample(ids,  length(ids),replace=TRUE)
  # Take all of the counties that each observation in our bootstrapped sample has 
  dataset <-  inner_join(tibble(GEOID = county_sample), senators_appropriations, by = "GEOID")
  
  # Then, run the zero-inflated model with this sample
  zim_2 <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                     core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                     freshman + pct_urban_scaled + median_household_income_scaled + factor(year) + capital,
                   zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                   data = dataset, family =  gaussian)
  
  summary_zim <- summary(zim_2)
  
  results_cond_coefs[[k]]  <- summary_zim$coefficients$cond[,1]
  results_cond_ses[[k]]  <- summary_zim$coefficients$cond[,2]
  results_zi_coefs[[k]]  <- summary_zim$coefficients$zi[,1]
  results_zi_ses[[k]]  <- summary_zim$coefficients$zi[,2]
  
  if (k %% 100 == 0) cat(paste(k, "out of", nboot, "samples stored..\n"))
  
}

# Turn the bootstrap results into our CI estimates
# Omitting models with NA SEs
to.keep <- which(apply(is.na(do.call(rbind.data.frame,results_cond_ses)),1,sum) == 0)

# SD of coefficients to create standard errors for bootstrap
cond_vec_se_alt <- apply(do.call(rbind.data.frame,results_cond_coefs)[to.keep,], 2, sd)
zi_vec_se_alt <- apply(do.call(rbind.data.frame,results_zi_coefs)[to.keep,], 2, sd)

# Variable names for our bootstrapped SEs
names(cond_vec_se_alt) <- c("intercept","log_countypop_scaled", "other_sen_requested", "other_sen_sameparty", "on_appropriations",
                            "female", "dem", "core_county","swing_county","seniority_scaled",
                            "party_leader", "meddist_scaled", "freshman", "pct_urban_scaled", 
                            "median_household_income_scaled", "factor(year)2022","capital","other_sen_requested:other_sen_sameparty")
names(zi_vec_se_alt) <- c("intercept","dem", "on_appropriations", "meddist_scaled", "dem:meddist_scaled")

# Rounding
cond_vec_se_alt <- round(cond_vec_se_alt,3)
zi_vec_se_alt <- round(zi_vec_se_alt,3)

# Running the zero inflated model outside of the loop to get coefficients
zim_2_boot2 <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary
summary_zim_boot2 <- summary(zim_2_boot2)

# Save coefficients from the model
cond_boot2 <- summary_zim_boot2$coefficients$cond
zi_boot2 <- summary_zim_boot2$coefficients$zi
coefs_boot2 <- c(cond_boot2[,1])
coefs2_boot2 <- c(zi_boot2[,1])

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model and the bootstrap above
basic_reg_boot2 <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + factor(year) + capital,
                data = senators_appropriations)

basic_reg2_boot2 <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# se takes the bootstrapped SEs
stargazer(basic_reg2_boot2, basic_reg_boot2,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from DW-NOMINATE Median",
                              "Distance from DW-NOMINATE Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "Fiscal Year 2023","County Includes Capital City", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. Standard errors come from 20000 bootstraps conducted at the county level. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_boot_county_level",
          digits=3,
          coef=list(coefs2_boot2, coefs_boot2), 
          se=list(zi_vec_se_alt, cond_vec_se_alt),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level with Bootstrapped Standard Errors (County-Level Bootstrap)",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))


##############
# TABLE C.13 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

zim_percapita <- glmmTMB(logpercap ~ log_countypop_scaled  + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_percapita <- summary(zim_percapita)
ses_percapita <-standard_error(zim_percapita)

# Save coefficients from the model
cond_percapita <- summary_zim_percapita$coefficients$cond
zi_percapita <- summary_zim_percapita$coefficients$zi
coefs_percapita <- c(cond_percapita[,1])
coefs2_percapita <- c(zi_percapita[,1])

# Save standard errors from the model
ses_cond_percapita <- c(ses_percapita$SE[ses_percapita$Component=="conditional"])
names(ses_cond_percapita) <- ses_percapita$Parameter[ses_percapita$Component=="conditional"]
ses_zi_percapita <- c(ses_percapita$SE[ses_percapita$Component=="zero_inflated"])
names(ses_zi_percapita) <- ses_percapita$Parameter[ses_percapita$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_percap <- lm(logpercap ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled  +
                         other_sen_requested*other_sen_sameparty +
                         female + core_county + swing_county + seniority_scaled + party_leader + 
                         freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                       data = senators_appropriations)

basic_reg2_percap <- lm(logpercap~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                        data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_percap, basic_reg_percap,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(Per Capita County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_percap",
          digits=3,
          coef=list(coefs2_percapita, coefs_percapita), 
          se=list(ses_zi_percapita, ses_cond_percapita),
          digits.extra = 0,
          title="Predictors of Per Capita Spending Requests at the County Level",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.14 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_2_divided_county <- glmmTMB(log_sum_by_county ~ log_countypop_scaled  + other_sen_requested*other_sen_sameparty + on_appropriations +  female + dem +
                                  core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                                  freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                                zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                                data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_divided_county <- summary(zim_2_divided_county)
ses_divided_county <-standard_error(zim_2_divided_county)

# Save coefficients from the model
cond_divided_county <- summary_zim_divided_county$coefficients$cond
zi_divided_county <- summary_zim_divided_county$coefficients$zi
coefs_divided_county <- c(cond_divided_county[,1])
coefs2_divided_county <- c(zi_divided_county[,1])

# Save standard errors from the model
ses_cond_divided_county <- c(ses_divided_county$SE[ses_divided_county$Component=="conditional"])
names(ses_cond_divided_county) <- ses_divided_county$Parameter[ses_divided_county$Component=="conditional"]
ses_zi_divided_county <- c(ses_divided_county$SE[ses_divided_county$Component=="zero_inflated"])
names(ses_zi_divided_county) <- ses_divided_county$Parameter[ses_divided_county$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_divided_county <- lm(log_sum_by_county ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled  +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                data = senators_appropriations)

basic_reg2_divided_county <- lm(log_sum_by_county ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_divided_county, basic_reg_divided_county,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_divided",
          digits=3,
          coef=list(coefs2_divided_county, coefs_divided_county), 
          se=list(ses_zi_divided_county, ses_cond_divided_county),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level, Total Appropriations Divided by the Number of Counties",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.15 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_2_continuous_safety <- glmmTMB(log_sum ~ log_countypop_scaled  + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   same_partyvote_scaled + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_continuous_safety <- summary(zim_2_continuous_safety)
ses_continuous_safety <-standard_error(zim_2_continuous_safety)

# Save coefficients from the model
cond_continuous_safety <- summary_zim_continuous_safety$coefficients$cond
zi_continuous_safety <- summary_zim_continuous_safety$coefficients$zi
coefs_continuous_safety <- c(cond_continuous_safety[,1])
coefs2_continuous_safety <- c(zi_continuous_safety[,1])

# Save standard errors from the model
ses_cond_continuous_safety <- c(ses_continuous_safety$SE[ses_continuous_safety$Component=="conditional"])
names(ses_cond_continuous_safety) <- ses_continuous_safety$Parameter[ses_continuous_safety$Component=="conditional"]
ses_zi_continuous_safety <- c(ses_continuous_safety$SE[ses_continuous_safety$Component=="zero_inflated"])
names(ses_zi_continuous_safety) <- ses_continuous_safety$Parameter[ses_continuous_safety$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_continuous_safety <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + same_partyvote_scaled + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                data = senators_appropriations)

basic_reg2_continuous_safety <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_continuous_safety, basic_reg_continuous_safety,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Previous County Voteshare","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_altvote",
          digits=3,
          coef=list(coefs2_continuous_safety, coefs_continuous_safety), 
          se=list(ses_zi_continuous_safety, ses_cond_continuous_safety),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Using a Continuous Measure of Electoral Safety",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.16 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_2_alt_cut <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   core_county2 + swing_county2 + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_alt_cut <- summary(zim_2_alt_cut)
ses_alt_cut <-standard_error(zim_2_alt_cut)

# Save coefficients from the model
cond_alt_cut <- summary_zim_alt_cut$coefficients$cond
zi_alt_cut <- summary_zim_alt_cut$coefficients$zi
coefs_alt_cut <- c(cond_alt_cut[,1])
coefs2_alt_cut <- c(zi_alt_cut[,1])

# Save standard errors from the model
ses_cond_alt_cut <- c(ses_alt_cut$SE[ses_alt_cut$Component=="conditional"])
names(ses_cond_alt_cut) <- ses_alt_cut$Parameter[ses_alt_cut$Component=="conditional"]
ses_zi_alt_cut <- c(ses_alt_cut$SE[ses_alt_cut$Component=="zero_inflated"])
names(ses_zi_alt_cut) <- ses_alt_cut$Parameter[ses_alt_cut$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_alt_cut <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county2 + swing_county2 + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + median_household_income_scaled + capital + factor(year),
                data = senators_appropriations)

basic_reg2_alt_cut <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_alt_cut, basic_reg_alt_cut,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County (57.5)", "Swing County (47.5)","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Median Household Income",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_altcoreswing",
          digits=3,
          coef=list(coefs2_alt_cut, coefs_alt_cut), 
          se=list(ses_zi_alt_cut, ses_cond_alt_cut),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Using Alternative Cutpoints for Core and Swing County",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.19 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_2_lowerthird <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + tercile_income + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family = gaussian)
# Model summary and standard errors
summary_zim_lowerthird <- summary(zim_2_lowerthird)
ses_lowerthird <- standard_error(zim_2_lowerthird)

# Save coefficients from the model
cond_lowerthird <- summary_zim_lowerthird$coefficients$cond
zi_lowerthird <- summary_zim_lowerthird$coefficients$zi
coefs_lowerthird <- c(cond_lowerthird[,1])
coefs2_lowerthird <- c(zi_lowerthird[,1])

# Save standard errors from the model
ses_cond_lowerthird <- c(ses_lowerthird$SE[ses_lowerthird$Component=="conditional"])
names(ses_cond_lowerthird) <- ses_lowerthird$Parameter[ses_lowerthird$Component=="conditional"]
ses_zi_lowerthird <- c(ses_lowerthird$SE[ses_lowerthird$Component=="zero_inflated"])
names(ses_zi_lowerthird) <- ses_lowerthird$Parameter[ses_lowerthird$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_lowerthird <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + tercile_income + capital + factor(year),
                data = senators_appropriations)

basic_reg2_lowerthird <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_lowerthird, basic_reg_lowerthird,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "Tercile Measure of Median Household Income (Low to High)",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_altincome",
          digits=3,
          coef=list(coefs2_lowerthird, coefs_lowerthird), 
          se=list(ses_zi_lowerthird, ses_cond_lowerthird),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Replacing Median Income with Tercile Measure of Median Household Income (Low to High)",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))

##############
# TABLE C.22 #
##############
# Remove all objects
rm(list=ls())
# Loading the data
load("county_level_request_data.RData")

# Running the zero inflated model
zim_2_poverty <- glmmTMB(log_sum ~ log_countypop_scaled + other_sen_requested*other_sen_sameparty + on_appropriations + female + dem +
                   core_county + swing_county + seniority_scaled + party_leader + meddist_scaled + 
                   freshman + pct_urban_scaled + per_poverty_scaled + capital + factor(year),
                 zi = ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations, family =  gaussian)
# Model summary and standard errors
summary_zim_poverty <- summary(zim_2_poverty)
ses_poverty <-standard_error(zim_2_poverty)

# Save coefficients from the model
cond_poverty <- summary_zim_poverty$coefficients$cond
zi_poverty <- summary_zim_poverty$coefficients$zi
coefs_poverty <- c(cond_poverty[,1])
coefs2_poverty <- c(zi_poverty[,1])

# Save standard errors from the model
ses_cond_poverty <- c(ses_poverty$SE[ses_poverty$Component=="conditional"])
names(ses_cond_poverty) <- ses_poverty$Parameter[ses_poverty$Component=="conditional"]
ses_zi_poverty <- c(ses_poverty$SE[ses_poverty$Component=="zero_inflated"])
names(ses_zi_poverty) <- ses_poverty$Parameter[ses_poverty$Component=="zero_inflated"]

# Running linear models to feed into stargazer, but we will replace the coefficients and SEs with those from the zero-inflated model above
basic_reg_poverty <- lm(log_sum ~ dem + on_appropriations + meddist_scaled + log_countypop_scaled +
                  other_sen_requested*other_sen_sameparty +
                  female + core_county + swing_county + seniority_scaled + party_leader + 
                  freshman + pct_urban_scaled + per_poverty_scaled + capital + factor(year),
                data = senators_appropriations)

basic_reg2_poverty <- lm(log_sum ~  dem + on_appropriations + meddist_scaled + dem*meddist_scaled,
                 data = senators_appropriations[complete.cases(senators_appropriations[,c("dem","on_appropriations","meddist_scaled","log_countypop_scaled",
                                                                                           "other_sen_requested","other_sen_sameparty",
                                                                                           "female","core_county","swing_county","seniority_scaled","party_leader", 
                                                                                           "freshman","pct_urban_scaled","per_poverty_scaled","capital","year")]),])

# Save the regression results as a table
# coef and se take the coefficients and standard errors from the zero-inflated model
stargazer(basic_reg2_poverty, basic_reg_poverty,
          omit=c("Constant"),
          notes.append = FALSE,notes.label = "",
          report="vc*s",star.char=c("*","**"),star.cutoffs = c(0.10,0.05),no.space = TRUE,
          font.size = "footnotesize", model.numbers = FALSE,
          column.labels=c("First Stage","Second Stage"),
          dep.var.labels = "Log(County Appropriations Requests + 1)",column.sep.width="0pt",
          covariate.labels=c( "Democrat (Majority Party Member)", "Member of Appropriations Committee","Distance from Floor Median",
                              "Distance from Floor Median x Democrat", "Log County Population",
                              "Other Senator Requested Funds to County","Other Senator is Same Party",
                              "Senator is a Woman","Core County", "Swing County","Seniority", "Party Leader", 
                              "Freshman Senator", "County Percent Urban Population", "County Percent Below the Poverty Line",
                              "County Includes Capital City", "Fiscal Year 2023", "Other Senator Requested * Other Senator Same Party"),
          notes="\\parbox[t]{.9\\textwidth}{\\footnotesize \\textit{Note}: The table shows the results from a zero-inflated gaussian regression of a
                        senator's county-level appropriation request behavior. The first stage models whether a senator makes a request at all or not and the
                        second stage models the logged total amount of funding a senator requests. $^{**}p<0.05$, $^*p<0.10$. All continuous variables are standardized for ease of comparison.}",
          label="tab1_county_poverty",
          digits=3,
          coef=list(coefs2_poverty, coefs_poverty), 
          se=list(ses_zi_poverty, ses_cond_poverty),
          digits.extra = 0,
          title="Predictors of Spending Requests at the County Level Replacing Median Income with Percent Below the Poverty Line",
          omit.stat = c("ll","rsq","adj.rsq","ser","f"))
